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I.  Summary 


A  study  of  optimal  aerodynamics  and  propulsive  control  at  supercircular  speeds  has  been 
Initiated.  The  objective  of  the  research  is  to  develop  methods  for  determining  optimal  guidance  ana 
control  of  earth  launch  kinetic  energy  weapons  designed  to  intercept  intercontinental  ballistic 
missiles  early  in  their  ascent  trajectory.  Optimal  control  techniques  are  being  used  to  obtain  multi¬ 
stage  trajectories  based  on  minimizing  the  mass  ratio.  Study  parameters  include  time  of  flight, 
down  range  intercept  distance  and  intercept  altitude.  Re-entry/skip  trajectories  are  considered. 
Innovative  means  of  attitude  control  of  the  final  stage  which  intercepts  the  target  are  being  studied. 

In  addition  an  investigation  of  the  control  of  aerodynamic  forces  on  hypersonic  vehicles  by 
boundary  layer  injection  has  been  started.  The  goal  of  this  work  is  to  determine  optimal  patterns  of 
injection  of  a  gas  into  a  boundary  layer  on  a  hypersonic  vehicle,  to  generate  desired  aerodynamic 
forces.  Two  directions  of  approach  are  being  studied.  In  the  first,  analytical  means  are  being  used 
to  study  the  effects  of  blowing  on  simple  flow  problems  in  the  various  flow  regimes;  a 
combination  of  asymptotic  and  numerical  methods  are  in  use.  In  the  second,  numerical  methods 
are  being  used,  with  particular  emphasis  on  obtaining  efficient  codes  which  result  in  the 
computation  of  crisp  shock  waves  and  which  can  handle  blowing  in  the  boundary  layer.  Both 
distributed  and  strip  blowing  are  under  consideration. 

Preliminary  results  are  presented. 

II.  Research  Objectives 

The  contract  for  this  research  work  has  two  phases,  one  dealing  with  optimal  trajectories  at 
superorbital  velocities  and  the  other  with  aerodynamic  force  calculations.  Hence  the  work 
statement  for  the  year's  work  has  two  separate  parts  as  follows,  where  paragraph  numbering  from 
the  contract  is  repeated  here. 
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t.  Formulate  the  automatic  control  problem  description  of  the  dynamics  of  ' 
transatmospheric  vehicles  for  ballistic  missile  defense.  Appropriate  control  laws  shall  be  L  ‘ 
determined  and  techniques  for  vehicle  optimization  shall  be  developed.  ----- 
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2.  Examine  schemes  for  the  generation  of  aerodynamic  control  forces  on  hypervelocity 
vehicles  by  the  injection  of  gas  into  the  boundary  layer.  Plan  analytical  and  numerical 
investigation  of  both  slot  and  surface  injection  from  hypersonic  slender  bodies  and  select 
techniques  for  distinguishing  among  dominant  fluid  mechanical  phenomena.  Minimize  the 
effects  of  numerical  difference  schemes  upon  shock  wave  thickness  and  location. 


III.  Status  of  the  Research  Effort 

This  work  began  on  August  15,  1986.  The  contract  was  to  be  for  37  months,  divided  into 
three  periods,  the  first  being  13  months  and  the  last  two  being  12  months  each.  At  the  end  of  the 
first  period,  pending  demonstration  of  satisfactory  work,  the  first  option  was  to  be  exercised  in 
order  that  work  during  the  second  period  could  proceed.  A  technical  briefing  was  given  to  Col.  D. 
Finkleman  (Air  Force),  Col.  W.  Merritt  (Army)  and  Major  S.  Rezwick  (Air  Force)  in  early  August 
of  1986;  it  was  indicated  that  they  were  well  satisfied  with  the  technical  effort  and  results. 
However,  during  the  time  that  option  one  could  have  been  exercised,  the  contract  was  being 
transferred  to  Army  supervision  and  apparently  the  group  transferring  the  contract  to  the  Army  was 
not  aware  that  further  action  was  required.  As  a  result  this  contract  is  terminated  after  13  months 
and  a  new  contract  for  the  remaining  two  options  is  being  renegotiated  with  the  Army  Strategic 
Defense  Command.  Thus,  this  final  report  covers  one  year  of  a  proposed  three  year  effort.  The 
following  status  report  is  in  two  pans,  with  the  numbering  corresponding  to  that  of  the  phases  in 
the  work  statement. 

1.  Optimal  Aerodynamic  and  Propulsive  Control  of  Intercept  Trajectories  at  Supercircular 
Speeds. 

Roben  M.  Howe,  Principal  Investigator 
Nguyen  X.  Vinh  Co-Investigator 
Elmer  G.  Gilbert  Co- Investigator 
Ping  Lu,  Graduate  Student 
Aharon  Karou,  Visiting  Scholar 
Gi!  Shorr,  Visiting  Scholar 

1,1  Introduction 

The  utilization  of  satellite-based  kinetic -energy  weapons  for  intercepting  intercontinental 
ballistic  missiles  during  their  ascent  trajectory  constitutes  one  of  the  major  efforts  in  the  current 
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U.S.  Strategic  Defense  Initiative  program.  In  this  research  we  have  considered  an  alternative  to 
space-based  weapons,  namely,  the  use  of  land  or  sea-based  kinetic  energy  interceptors.  When 
such  interceptors  are  based  on  the  U.S.  continent,  ranges  of  up  to  5000  miles  are  needed.  If  the 
interceptors  are  sea  launched,  for  example  from  submarines,  the  required  interceptor  range  can  be 
reduced  to  as  little  as  2000  miles  or  less.  In  either  case  the  interceptor  must  accelerate  to  very  high 
speeds,  typically  twice  circular  orbit  velocity,  in  order  to  reach  a  target  at  these  distances  within 
several  minutes. 

In  order  to  achieve  these  velocities,  a  rocket  powered  interceptor  requires  an  overall  mass 
ratio  of  several  thousand  to  one.  Even  so,  the  interceptor  takeoff  weight  can  be  held  within 
reasonable  bounds  by  the  use  of  miniature  guided  warheads.  It  is  clear  that  such  large  mass  ratios 
can  only  be  achieved  by  using  between  4  and  7  rocket  boost  stages.  In  this  research  we  have 
considered  the  problem  of  optimizing  a  number  of  stage  parameters  and  the  angle  of  attack  control 
in  order  to  minimize  the  overall  required  mass  ratio.  The  time  history  of  angle  of  attack  was 
parameterized  by  letting  the  angle  of  attack  for  each  thrusting  stage  vary  linearly  with  time. 
Additional  parameters  to  be  optimized  included  the  thrust  bum  times  for  each  stage  and  the  coast 
times  between  stage  bums.  Although  rockets  of  this  type  would  normally  be  launched  straight  up, 
we  have  also  included  the  initial  launch  flight-path  angle  as  a  parameter.  This  was  based  on  the 
assumption  that  a  fairly  rapid  pitchover  maneuver  can  be  performed  immediately  after  launch  to 
achieve  the  desired  initial  flight  path  angle. 

Once  the  trajectory  optimization  problem  has  been  converted  to  a  parameter  optimization 
problem,  a  quasi-Newton  minimization  algorithm  is  utilized  with  penalty  functions  to  implement 
the  terminal  constraints  required  for  intercept,  as  well  as  any  additional  constraints.  The  neces  'try 
gradients  of  the  cost  function  with  respect  to  the  parameters  are  computed  numerically  using  tn  w 
differences.  A  complete  multi-stage  trajectory  must  be  run  using  numerical  integration  for  ev  .  v 
evaluation  of  the  cost  function.  For  each  set  of  N  parameters  this  means  that  N+l  con.p  * 
trajectories  must  be  computed,  one  to  evaluate  the  cost  function  for  the  given  parameter  settings, 
and  N  to  evaluate  the  gradient  of  the  cost  function  with  respect  to  each  of  the  N  para' 
the  trajectory  optimization  problems  we  have  considered,  between  10  and  20  parameu  -*■ -n 

used.  Thus  between  1 1  and  21  complete  multi-stage  trajectories  must  be  computed  nuirtatd,iy 
a  single  iteration  of  the  optimization  procedure.  Also,  each  iteration  utilizes  a  line  search  which 
requires  several  additional  trajectory  computations.  It  has  been  found  that  between  200  and  400 
iterations  are  needed  until  satisfactory  convergence  to  an  optimal  solution  is  achieved.  It  follows 
that  several  thousand  trajectories  must  be  computed  for  each  optimization  case  considered.  In 
order  to  reduce  the  overall  computational  time,  the  multiple-stage  trajectories  were  calculated  using 


3 


an  AD  100  multiprocessor  computer.  A  VAX  host  computer  was  used  to  implement  the 
optimization  algorithm  and  input  new  parameters  to  the  AD  100  for  each  successive  trajectory  run. 
In  this  way  the  total  time  required  for  an  overall  trajectory  optimization  was  reduced  from  several 
hours  to  several  minutes. 

In  the  next  section  scalar  equations  of  motion  for  the  interceptor  trajectory  are  presented. 
Subsequent  sections  discuss  the  optimization  problem,  its  numerical  solution,  and  results  for  a 
typical  4-stage  trajectory  optimization. 

1.2  Equations  of  Motion 

In  order  to  reduce  the  trajectory  computation  times,  only  the  two  dimensional  equations  of 
motion  of  the  interceptor  were  considered.  As  a  further  simplification,  the  effect  of  earth's  rotation 
was  neglected.  It  was  felt  that  neither  of  these  simplifications  altered  the  overall  nature  of  the 
optimization  problem.  More  comprehensive  equations  of  motion  can  always  be  incorporated  into 
the  optimization  problem  at  a  later  time  to  improve  the  accuracy  for  those  cases  which  may  turn  out 
to  be  of  practical  interest.  Thus  the  interceptor  rocket  was  considered  to  be  a  point  mass  m  moving 
in  a  plane  which  contains  the  center  of  a  spherical,  non-rotating  earth  with  inverse  square 
gravitational  field.  A  convenient  and  efficient  method  for  representing  the  velocity  state  of  the 
interceptor  is  in  terms  of  the  total  velocity  and  the  flight  path  angle  with  repect  to  the  local 
horizontal  [1].  This  leads  to  the  following  two-dimensional  equations  of  motion: 
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where 

A,  =  rocket  length  d  -  rocket  diameter 

The  atmospheric  density,  p,  in  Eq.  (1.6)  is  a  function  of  the  rocket  altitude  h,  which  is 
given  by  the  equation 
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h  -  (R  -  1)  ro 


(1.9) 


Three  different  methods  for  approximating  the  variation  of  p  with  h  were  considered.  The  first 
utilizes  table  lookup  with  linear  interpolation  based  on  a  tabulated  standard  atmosphere.  This 
method  introduces  into  the  simulation  a  function  with  discontinuous  slopes,  which  has  the  potential 
of  causing  numerical  difficulties  in  the  descent  algorithm  (see  Section  1 .4).  The  second  method 
uses  the  following  exponential  approximation: 


P  = 


pO  e  '  Ph 


(1.10) 


where  1/p  =  23,500  ft.  and  p0  is  the  atmospheric  density  at  sealevel.  The  third  method  uses  the 
following  analytic  approximation  [3]: 


p/pO  =  [Ao  +  Ajh  +  A2h2  +  •  •  •  •  +  Ajihl  1]  ’  4 


(1.11 


with 

Ao  =  1.0000000000 
A3  =  0.5497466428  x  10*3 
A6  = -0.2291755793  x  10*7 
A9  =  0.1010575266  x  10‘13 


Aj  =0.3393495800  x  10'1 
a4  = -0.3228358326  x  10*4 
A7  =0.2902146443  x  10*9 
A  io=  -0.2482089627  x  10* 16 


A2  = -0.3433553057  x  10' 2 
a5  =  0.1106617734  x  10' 5 
A8  =  -0.2230070938  x  10*11 
An  =0.2548769715  x  10‘19 


Eq.  (1.11)  is  valid  from  0  to  200  km.  altitude  and  is  accurate  to  better  than  5  percent  for  altitudes 
up  to  70  km.  It  is  a  considerably  better  approximation  than  the  exponential  model  in  Eq.  (MO) 
and  takes  no  longer  for  computer  execution  than  the  table  lookup  scheme.  Thus  it  represented  our 
preferred  approach. 


6 


Dimensionless  versions  of  Eqs.  (1.1)  through  (1.9)  and  Eq.  (1.11)  were  used  for  the 
trajectory  simulation  of  each  inteceptor  stage  while  thrusting  or  coasting  (T  =  0).  Tne  state 
variables  r,  0,  V  and  y  are  continuous  from  one  stage  to  the  next.  On  the  other  hand,  the  state 
variable  m,  which  represents  the  mass  of  each  stage,  undergoes  a  discontinuous  decrease  each  time 
a  stage  bums  out  and  is  dropped.  The  mass  state  m  can  be  made  continuous  by  defining  it  as  the 
mass  of  rocket  fuel  divided  by  the  initial  takeoff  mass.  The  totai  instantaneous  mass  for  each  stage 
is  then  the  sum  of  the  fuel  mass  and  the  inert  mass.  The  bum  time  for  each  stage,  as  well  as  the 
coast  time  per  stage,  constitute  parameters  to  be  optimized.  For  a  given  stage  mass  ratio,  the  stage 
bum  time  determines  the  thrust  level  for  the  stage.  Thus  the  parameters  representing  stage  bum 
times  are  equivalent  to  parameters  representing  stage  thrust  levels.  As  noted  previously,  the  time 
history  of  angle  of  attack  was  parameterized  by  letting  it  vary  linearly  for  each  stage.  Thus  an 
additional  parameter  per  stage  is  the  time-rate-of-change  of  angle  of  attack. 

1.3.  Constraints  and  Optimization  Problem  Formulation 

To  complete  the  statement  >r  the  optimization  problem  it  was  necessary  to  introduce 
certain  constraints  on  the  trajectory  and  its  defining  parameters.  The  most  obvious  of  these 
were  the  ones  required  for  target  interception.  After  final  stage  cutoff,  the  payload  must 
coast  in  such  a  manner  that  at  the  specified  intercept  time,  tf,  the  position  variables,  r(tf)  and 
0(tf),  match  those  of  the  target. 

There  are  several  ways  to  implement  the  intercept  conditions.  One  way  is  to  impose  them 
at  the  time  of  final-stage  cutoff,  tc.  If  r(tc),  0(tc),  r(tf),  and  0(tf)  are  known  and  the  coast  is  an  out- 
of-the-atmosphere  conic  transfer,  the  values  yftc)  and  V(tc)  required  for  intercept  are  determined  by 
solving  an  appropriate  Lambert  problem.  This  can  be  done  numerically  by  the  methods  given  in 
[4],  The  errors  between  the  required  and  actual  values  of  7^)  and  VltO  must  be  zero  for  intercept 
and  these  conditions  give  two  equality  constraints  on  the  trajectory  at  t  =  tc. 

It  was  found  that  unconstrained  optimal  conic  coast  trajectories  e-enter  the  atmosphere  or 
even  pass  through  the  earth's  surface.  To  assure  that  this  doesn't  happen,  the  minimum  altitude 
during  the  coast  was  computed  (as  part  of  the  solution  of  the  L  ambert  problem)  and  specified  to  lie 
outside  the  sensible  atmosphere.  The  error  between  the  computed  minimum  altitude  and  its 
specified  value  constituted,  when  set  to  zero,  a  third  equality  constraint. 

Another  way  to  implement  the  intercept  condition  is  to  integrate  numerically  the  equations 
of  motion  of  the  final  coast  trajectory,  just  as  the  equations  for  the  previous  stages  were  integrated. 
In  this  case  the  errors  between  the  interceptor  coordinates  rt  tf),  0(  tf)  and  the  target  coordinates  at  the 
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final  time  tf  are  set  equal  to  zero  to  form  two  equality  constraints.  As  before,  a  minimum  allowable 
altitude  outside  the  sensible  atmosphere  can  be  chosen  for  the  final  coast  trajectory.  This  minimum 
altitude  is  then  the  third  constraint.  Alternatively,  the  effect  of  the  atmosphere  itself,  which  is 
included  in  the  trajectory  equations,  will  ensure  that  a  realistic  altitude  constraint  is  enforced  when 
the  final  coast  trajectory  equations  arc  integrated.  Indeed,  this  is  the  required  approach  in  the  aero- 
assisted  case,  where  the  final  coast  trajectory  actually  utilizes  aerodynamic  down  lift  to  hold  the 
trajectory  at  near  constant  altitude  in  the  presence  of  supercircular  velocity.  In  Section  1.5  an 
example  of  this  type  is  presented. 

In  addition  to  the  equality  constraints  on  the  trajectory,  inequality  constraints  may  be 
imposed  on  some  of  the  problem  parameters.  For  example,  an  initial  stage  acceleration,  A«q\  may 
have  upper  and  lower  limits,  Amax  and  Amjn.  Direct  treatment  of  these  constraints  in  the 
algorithmic  process  was  avoided  by  imposing  them  implicitly  through  the  introduction  a  nonlinear 
transformation  of  a  new  parametric  variable.  In  the  case  just  mentioned  a  suitable  transformation  is 
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where  ay  is  the  new  unconstrained  parameter. 

Putting  all  of  the  above  details  together  in  a  compact  notation  yields  a  problem  of  the  form 
minimize  f(z)  subject  to  g,(z)  =  0  ,  i  =  1 ,2,3.  (1.13) 

Here:  f(z)  is  tire  ratio  of  the  launch  mass  to  the  payload  mass,  the  g,(z)  are  the  errors  in  meeting  the 
equality  constraints,  and  z  is  a  vector  whose  components  are  the  "unconstrained"  parameters 
arising  from  models  described  in  Section  1.2.  Given  z,  it  is  clear  that  f(z)  and  gi(z)  may  be 
evaluated  by  integrating  the  equations  of  motion,  computing  mite),  yU*.-),  V(tc),  rit^),  and  0(te ).  and 
solving  the  Lambert  problem  determined  by  the  target  specification.  When  the  alternative  method 
involving  numerical  integration  of  the  final  coast  trajectory  is  used  to  implement  the  intercept 
condition,  f(z)  was  still  evaluated  from  m(tc),  but  gi(z)  is  evaluated  from  rf tf)  and  0(tf). 

The  algorithm  for  the  numerical  solution  of  (1.13)  requires  the  gradient  (collection 
of  N  first  partial  derivatives)  of  functions  such  as  f(z).  While  it  is  possible  to  derive  formulas  for 
the  partials  (they  may  be  written  in  tenns  of  the  differential  equations  adjoint  to  the  linearized 
equations  of  motion),  they  are  exceedingly  complex  because  of  the  nonlinear  functions  appearing 
in  Eqs.  (1.1)  -  (1.9)  and  (111),  llus  was  the  reason  for  using  finite  differences  for  the  gradient 
calculations. 


8 


1.4  Numerical  Solution  of  the  Optimization  Problem 


A  variety  of  numerical  approaches  have  been  used  to  solve  trajectory  optimization  problems 
with  equality  constraints.  These  include:  (1)  gradient  projection  methods,  (2)  Newton-Raphson 
solution  of  the  necessary  conditions,  (3)  penalty  function  me 'hods.  Approach  (1)  gives  slow 
convergence;  approach  (2)  requires  second  derivatives  of  problem  data  and  good  initial  estimates; 
approach  (3)  leads  to  badly  conditioned  (though  unconstrained)  minimization  problems.  In  the 
research  reported  here  approach  (3)  has  been  used,  with  its  disadvantages  attacked  by  using  a 
rapidly  convergent  quasi-Newton  algorithm  [5,6]  and  an  augmented  Lagrangian  for  the  penalty 
function. 

In  the  usual  penalty  function  method  approximate  solutions  of  (1.13)  are  sought  by 
algorithmic  minimization  of 

3 

Fp(z)  =  f(z)  +  £  ^lg,(z)J2  (L14) 

i  =  l 


It  can  be  shown  [5,6]  that  the  equality  conditions  can  be  accurately  enforced  by  choosing  the 
penalty  coefficients,  >  0  .  to  be  sufficiently  large.  In  practice,  numerical  hazards,  such  as  slow 
convergence  or  convergence  'O  false  minima,  are  reduced  by  minimizing  Fp(z)  several  times  with 
successively  larger  values  for  the  penalty  coefficients.  At  each  successive  stage  the  starting  point 
for  the  minimization  algorithm  is  taken  as  the  solution  point  from  the  previous  stage.  Even  so,  the 
conditioning  of  Fp(z)  worsens  [5,6]  as  the  p,  increase  and  there  is  a  practical  limit  to  the  accuracy 
which  can  be  achieved.  Numerical  errors  in  the  minimization  algorithm  and  the  evaluation  of  FpU) 
and  its  gradient  eventually  become  the  limiting  factors. 

In  the  augmented  Lagrangian  approach  Fp(z)  is  replaced  by 


Fl(z) 


=  f (z)  + 


A  g  (z) 

I  °\ 


1  =  1 


t 


F,[g,L 


(1.15) 


where  the  Aj  are  real  numbers  approximating  the  Lagrange  multipliers  for  the  minimization  problem 
(1.13).  When  F"l(z)  is  minimized,  it  is  possible  to  satisfy  the  equality  constraints  accurately  with 
much  smaller  values  for  the  penalty  coefficients  [5,6).  The  resulting  improved  conditioning  of 
1:I  (z)  reduces  the  bad  effects  of  the  numerical  errors.  Again,  it  is  advantageous  to  solve  a 
sequence  of  minimization  problems  with  successive  initializations,  where  in  this  case  both  the  A., 
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and  the  |ii  are  adjusted.  The  adjustment  scheme  proposed  by  Powell  and  described  on  page  134 
of  [5]  was  used  with  good  results  in  the  research  reported  here.  It  eliminated  the  need  for  human 
intervention  in  selecting  the  parameters  Xj  and  |ij  and  improved  by  at  least  an  order  of  magnitude 
die  accuracy  of  equality  constraints  over  the  simple  penalty  function  approach. 

Because  of  die  highly  nonlinear  character  of  Fl(z)  and  its  poor  conditioning,  it  was 
necessary  to  pay  careful  attention  to  the  procedure  for  its  unconstrained  minimization.  Certainly, 
the  second  order  properties  of  Fl(z)  must  be  taken  into  account.  The  BFGS,  quasi-Newton 
program  of  Shanno  and  Phua  [7]  was  found  to  be  effective.  Its  line  search  was  modified  to 
compute  derivatives  along  the  st  rch  direction  directly  by  forward  differences.  This  eliminated  the 
need  for  gradient  evaluations  during  the  line  search. 

The  performance  of  a  quasi-Newton  algorithm  is  predicated  on  the  smoothness  of  the 
function  being  minimized.  The  function  should  be  at  least  twice  continuously  differentiable.  Here, 
the  representation  of  empirical  data  appearing  in  the  equations  of  motion  is  important.  See,  for 
example,  the  comments  in  Section  1.2  on  computing  the  atmospheric  density  b>  table  lookup  with 
linear  interpolation. 

Another  concern  is  the  forward  difference  computation  of  the  partial  derivatives  of  Fl(z). 
If  the  step  size  is  too  small,  the  effects  of  round  off  errors  are  exacerbated;  if  the  step  size  is  too 
large,  truncation  errors  become  appreciable.  If  the  machine  precision  is  e,  it  can  be  argued  [6]  that 
the  step  size  should  be  the  order  of  e1/2.  In  our  case  the  machine  precision  e  @  10' 12  was 
determined  by  the  AD  100,  which  has  a  40  bit  significant!.  It  was  found  that  the  step  size  iO"6  was 
reasonable,  but  that  experimentation  with  the  step  size  could  improve  the  accuracy  noticeably. 
Even  better  accuracy  was  obtained  by  using  central  differences  for  the  partial  derivatives.  It 
reduced  the  effects  of  round-off  errors,  since  in  this  case  a  larger  step  size  could  be  used  for  the 
same  truncation  error  16],  Of  course,  2N  +  1  evaluations  of  F'l(z)  are  then  needed  to  obtain  Fl(z) 
and  its  gradient. 

In  general  the  trajectory  state  equations  were  solved  numerically  using  fourth-order  Runge- 
Kutta  integration  with  a  fixed  step  size.  A  single-pass  predictor  method  with  a  Gear  RK-4  startup 
algorithm  was  found  to  give  better  peiformance  in  some  cases  [8]  and  represents  an  area  for  further 
research. 


Two  example  optimal  trajectory  solutions  are  considered  in  this  section.  0(tf)  -  0 q,  the  total 
change  in  polar  angle  of  each  trajectory,  is  set  equal  to  30  degrees.  This  corresponds  to  an 
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interceptor  range  of  2094  miles.  The  total  time  of  interceptor  flight  was  prescribed  to  be  300 
seconds  and  the  target  was  intercepted  at  an  altitude  of  200  km.  Stage  takeoff  weight  was  equal  to 
10,000  kg  and  physical  parameters  were  assumed  to  be  similar  to  those  of  a  Minuteman  III.  The 
specific  impulse  Isp  of  the  rocket  engines  was  equa*  to  300  seconds,  and  the  r  atio  of  initial  fuel 
weight  to  inert  weight  for  the  rocket  motor  of  each  stage  was  assumed  to  be  10.  The  intercept 
condition  was  implemented  by  two  equality  constraints  that  required  the  interceptor  coordinates  at 
the  final  time  tf  to  match  the  target  coordinates.  Two  cases  were  considered;  aeroassisted  and  non- 
aeroassisted.  In  the  aeroassisted  case  the  final  coasting  stage  reenters  the  earth's  atmosphere  with  a 
controlled  angle  of  attack  which  could  generate  downward  lift.  In  the  non-aeroassisted  case  the 
angle  of  attack  of  the  final  coasting  stage  w'as  maintained  at  zero,  which  meant  that  aerodynamic 
forces  could  only  cause  drag.  In  each  case  there  were  4  booster  stages,  with  equal  mass  ratios 
prescribed  for  each  stage.  The  common  mass  ratio  was  therefore  a  parameter.  As  noted  earlier, 
the  performance  objective  was  to  maximize  the  final  payload  weight.  The  results  are  summarized 
in  Table  1.1. 


TABLE  1.1.  Optimal  Trajectory  Data  for  4-stage  Interceptor 


Final  optimal  payload  weight 

Aeroassisted 

9.73  kg 

Non-aeroassisted 
5.63  kg 

Mass  ratio  per  stage 

3.98 

4.33 

Minimum  altitude  in  final  coast  to  target 

93.82  km 

96.82  km 

Burnout  speed  (Vc  =  circular  orbit  speed) 

1.766  Vc 

1.775  Vc 

Burnout  altitude 

98.24  km 

155.04  km 

Launch  flight  path  angle,  y0 

82.8  deg 

60.6  deg 

Initial  thrust  acceleration,  stage  1 

13.05  g 

9.28  g 

Bum  time,  stage  1 

17.21  sec 

24.86  sec 

Coast  time  between  stages  1  and  2 

6.49  sec 

1.84  sec 

Initial  thrust  acceleration,  stage  2 

23.37  g 

12.93  g 

Bum  time,  stage  2 

9.61  sec 

17.84  sec 

Coast  time  between  stages  2  and  3 

5.35  sec 

2.41  sec 

Initial  thrust  acceleration,  stage  3 

31.40  g 

23.87  g 

Bum  time,  stage  3 

7.15  sec 

9.67  sec 

Coast  time  between  stages  3  and  4 

14.85  sec 

56.74  sec 

Initial  thrust  acceleration,  stage  4 

37.52  g 

31.67  g 

Bum  time,  stage  4 

5.99  sec 

7.29  sec 

Final  coast  time  to  target 

233.67  sec 

179.36  sec 

Total  flight  time 

300  sec 

300  sec 

Total  distance 

2094  miles 

2094  miles 

11 


From  the  table  it  is  clear  that  the  optimal  aeroassisted  trajectory  yields  a  significantly  larger 
payload,  9.73  kg  versus  5.63  kg  for  the  non-aeroassisted  case.  In  Figure  1.1  are  shown  time 


Figure  1.1.  Time  history  plots  for  the  aeroassisted  optimal  trajectory. 
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history  plots  of  velocity,  flight  path  angle,  dynamic  pressure  and  angle  of  attack  for  the 
acroassisted  optimal  trajectory.  Figure  1.2  shows  altitude  versus  polar  angle  plots  for  both  the 
aeroassisted  and  non-aeroassisted  trajectories.  Note  how  the  acroassisted  trajectory  is  able  to 
maintain  an  approximately  constant  altitude  of  some  95  km  over  more  titan  half  of  the  final  coast  to 
intercept.  Th  !:?  possible  despite  the  supcrcircular  speed  because  of  the  aerodynamic  down  lift, 
which  rcac  naximum  of  1.87  g  during  this  portion  of  the  trajectory.  The  maximum  total 
downward  s  ;  itti^n,  including  gravity,  is  therefore  2.87g.  Since  the  centrifugal  acceleration  in 
g  units  is  V  *•  ,  is  the  circular-orbit  velocity),  it  follows  that  constant  altitude  (i.e.,  a  circular 
orbit)  will,  >.n  c:  presence  of  1.87  g  of  down  lift,  be  maintained  at  a  velocity  of  (2.87)1/2  or 
1.694VC.  Rsfci  ice  to  the  plot  of  V  versus  time  in  Figure  1.1  shows  that  this  is  indeed  the  case. 
On  the  other  ha,  we  see  in  Figure  1.2  that  the  non-aeroassisted  trajectory  must  climb  to  over  150 
km  before  the  d  1  stage  bum  in  order  to  avoid  significant  penetration  of  the  atmosphere.  This 
results  in  a  co  table  more  expensive  trajectory  in  terms  of  required  mass  ratio. 


Figure  1 .2  Comparison  of  optimd  tractories  with  and  without  aeroassistance. 
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The  plot  of  angle  of  attack  a  versus  time  in  Figure  1.1  shows  clearly  the  linear  variation  of 
oc  for  each  stage.  Note  that  zero  angle  of  attack  was  assumed  during  the  coast  between  stage 
bums.  This  was  based  on  the  assumption  that  the  interceptor  is  aerodynamically  stable  during 
these  coasts  with  no  active  thrust  vector  control.  Means  of  active  attitude  control  would  be 
required  to  maintain  the  linear  angle  of  attack  variation  shown  for  the  final  coast  to  intercept.  One 
option  for  achieving  this  would  be  to  provide  thrust  vector  control  utilizing  a  low  level  of  thrust 
during  this  final  portion  of  the  trajectory.  Such  thrust  will  probably  be  required  in  any  event  to 
provide  course  corrections  to  intercept  the  target.  When  this  is  included  in  the  overall  trajectory 
optimization,  it  may  very  well  produce  a  non-aeroassisted  optimal  trajectory  which  is  competitive 
with  the  aeroassisted  case,  since  now  a  downward  component  of  thrust  can  be  used  to  teplace  the 
aerodynamic  down  lift.  It  should  be  noted  that  boundary  layer  injection  at  hypersonic  speeds, 
which  is  the  subject  of  Part  2  of  the  research  reported  here,  is  also  a  candidate  for  attitude  control 
while  the  interceptor  is  either  coasting  or  thrusting  in  the  atmosphere. 

1 .6  Optimal  Intercept  from  a  Space  Base 

The  final  stage  of  the  computed  trajectory  is  essentially  a  small  rocket  interceptor  launched 
from  space,  i.e.,  from  the  ascending  (n-l)th  stage.  Therefore,  an  analytic  study  of  the  optimal 
intercept  trajectory  has  been  made  for  the  case  of  space  launch  from  an  orbiting  platform.  The 
launch  base  is  assumed  to  be  on  a  Keplerian  orbit.  The  trajectory  of  the  target  is  arbitrary,  but  it  is 
assumed  to  be  well  determined  after  a  time  to  called  the  acquisition  time.  Again,  the  intercept  time 
tf-to  is  restricted  to  several  minutes. 

For  flight  in  a  vacuum  and  in  a  central  Newtonian  force  field  with  impulsive  thrust,  a 
complete  analytical  solution  has  been  obtained.  The  study  used  the  well  established  theory  of  the 
primer  vector  in  optimal  transfer.  For  a  specified  intercept  time  tf,  the  one-impulse  solution 
initiated  at  the  time  to  is  assumed.  The  associated  Lambert  problem  is  solved  and  the  magnitude  of 
the  required  impulse  is  computed.  The  resulting  initial  conditions  and  the  transversality  conditions 
for  optimality  are  sufficient  to  compute  the  primer  vector  which  governs  the  optimal  thrust  control. 
Then,  based  on  the  information  provided  by  the  primer  vector,  rules  have  been  established  to 
search  for  the  optimal  solution  if  the  assumed  initial  one-impulse  trajectory  is  rot  optimal.  It  has 
been  found  that  there  are  three  possible  optimal  trajectories. 

a.  One  impulse  trajectory  with  the  thrust  applied  immediately  at  the  acquisition  time  tq. 

b  One-impulse  with  an  initial  coasting  phase  until  the  optimal  time  ti>to  before  the 
application  of  the  impulse. 
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c.  Two-impulse  trajectory  with  the  first  inpulse  at  to  and  the  second  impulse  at  mid-course 
for  final  interception. 

The  theory  can  be  used  to  explain  the  necessary  coast  time  between  the  nth  stage  of  the 
previously  computed  optimal  trajectories  for  a  multi-stage  interceptor  launched  from  the  surface  of 
the  earth. 

A  numerical  example  has  been  used  to  test  the  theory.  The  interceptor  is  launched  from  a 
base  in  a  circular  orbit  at  600  km  altitude.  The  target  is  intercepted  at  an  altitude  of  250  km  with  a 
down  range  longitude  of  45  degrees  and  various  latitudes  ranging  from  0  to  60  degrees.  The 
intercept  times  range  from  2  minutes  to  15  minutes.  For  each  case  the  maximum  characteristic 
velocity  for  intercept  is  computed.  Depending  on  the  prescribed  intercept  time,  the  trajectory  may 
change  from  an  elliptic  trajectory  to  a  clearly  hyperbolic  trajectory. 

1.7  Guidance  and  Control  to  Intercept  of  the  Final  Stage 

To  provide  terminal  control  out  of  the  atmosphere  so  the  final  stage  can  home  on  the  target 
requires  some  type  of  guidance  law  and  trajectory  control.  The  trajectory  can  be  changed  by 
providing  acceleration  at  right  angles  to  the  interceptor  flight  path.  For  an  interceptor  which  utilizes 
a  fixed  rocket  engine  to  provide  acceleration  along  its  longitudinal  axis,  a  required  lateral 
acceleration  component  can  be  generated  by  the  appropriate  change  in  attitude  angle. 

To  minimize  weight  and  complexity,  an  attitude  control  method  which  uses  only  a  single 
small  control  jet  has  been  studied.  This  single  jet  provides  thrust  at  right  angles  to  the  longitudinal 
axis  of  an  axially-symmetric  terminal  stage.  The  stage  is  given  a  large  roll  rate  about  its  longitudinal 
axis,  along  which  the  main  rocket  engine  provides  continuous  thrust.  Attitude  control  is  achieved 
by  turning  ort  the  small  side  jet  for  a  fraction  of  each  revolution  and  at  the  right  time  during  the  roll 
cycle.  A  simple  incremental  change  in  attitude  of  the  spinning  stage  can  be  viewed  as  a  change  in 
direction  of  its  angular  momentum  vector  without  a  change  in  vector  magnitude.  Each  discrete 
change  in  attitude  requires  at  least  two  impulses  from  the  control  jet.  The  first  impulse  produces  an 
increment  in  angular  momentum  at  right  angles  to  the  spin  axis.  The  second  control  jet  impulse  is 
applied  half  a  roll  period  later  to  cancel  the  angular  momentum  increment  at  right  angles  to  the  spin 
axis,  which  has  now  been  changed  in  direction  by  the  attitude  angle  increment.  The  size  of  the 
incremental  change  in  attitude  angle  can  be  varied  by  changing  the  thrust  duration  of  the  control  jet. 
If  the  available  thrust  is  insufficient  to  change  the  attitude  by  the  desired  amount,  then  a  number  of 
cyclic  impulses  followed  by  the  same  number  of  cyclic  impulses  initiated  half  a  period  later  can  be 
utilized. 
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The  dynamics  of  the  above  control  scheme  has  been  studied  using  both  analytic  methods 
and  computer  simulation.  Minimum  resolution  in  attitude  angle  has  been  established  as  a  function 
of  system  parameters,  including  minimum  control -jet  thrusting  time.  Multiple  thrust  maneuvers  for 
large  attitude  changes  have  also  been  studied  in  the  same  manner.  Required  control  laws  for 
achieving  desired  attitude  angle  changes  have  been  developed. 

Once  a  means  of  attitude  control  is  available,  a  guidance  law  is  needed  to  specify  the 
required  trajectory  changes.  The  most  commonly  used  guidance  law  for  homing  missiles  is 
proportional  navigation.  Here  it  does  not  appear  to  be  a  practical  guidance  method  because  it 
results  in  excessive  terminal  accelerations  for  the  interceptor.  As  an  alternative,  a  simple  algorithm 
has  been  developed  to  estimate  the  future  positions  of  both  interceptor  and  target.  Based  on  target 
position,  velocity  and  acceleration,  and  the  known  equations  of  motion,  a  power  series  solution  for 
future  target  position  is  obtained.  In  the  same  way,  the  future  interceptor  position  is  computed. 
From  these  predicted  trajectories  the  required  interceptor  attitude  for  zero  miss  is  calculated  and 
utilized  to  implement  the  guidance  law.  Computer  simulations  have  been  used  to  test  the 
effectiveness  of  this  method  in  combination  with  the  single-axis  attitude  control  scheme  described 
above.  The  simulations  have  been  run  for  both  thrusting  and  coasting  targets.  Based  on 
reasonable  interceptor  parameters,  including  roll  rates,  control-thrust  levels,  and  minimum  control- 
thrust  times  it  has  been  determined  that  direct  hits  can  be  obtained.  However,  these  simulations 
have  not  yet  taken  into  account  the  effect  of  guidance  measurement  errors,  nor  have  they 
considered  the  effect  of  target  thrust  termination  during  terminal  guidance  of  the  interceptor. 

1.8  Conclusions  and  Subsequent  Research 

This  research  has  considered  techniques  for  computing  the  optimal  trajectories  of  earth- 
launched  interceptor  rockets  which  are  accelerated  to  over  twice  orbital  velocity.  The  objective  of 
the  optimization  is  to  maximize  the  payload  with  a  given  interceptor  takeoff  weight.  Terminal 
constraints  have  been  enforced  using  penalty  functions.  The  required  gradients  of  the  cost 
functions  have  been  computed  numerically  using  finite  differences.  Results  have  been  presented 
for  a  4-stage  example.  It  was  found  that  allowing  a  coast  period  between  booster  stages  and  using 
an  aeroassisted  trajectory  to  generate  down  lift  during  the  final  coast  segment  improves  the  optimal 
payload  significantly.  The  research  has  also  considered  attitude  control  of  a  final  interceptor  stage 
which  uses  a  fixed  rocket  engine  to  provide  thrust  along  its  longitudinal  axis,  about  which  the 
rocket  also  spins.  Attitude  changes  are  achieved  by  firing  a  single  small  thruster  at  right  angles  to 
the  spin  axis  at  the  appropriate  time  during  each  spin  cycle.  In  combination  with  a  guidance  law 
based  on  predicted  future  target  and  interceptor  trajectories,  it  has  been  found  through  computer 
simulation  that  direct  hits  can  be  obtained. 
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Future  optimal  trajectory  research  efforts  will  include  parametric  studies  of  the  effect  of 
different  prescribed  flight  times,  target  range  and  altitude,  payload  mass,  and  refined  aerodynamic 
and  propulsion  models.  The  effect  of  aerodynamic  heating  on  the  ascent  trajectory  and  during  the 
aeroassisted  terminal  trajectory  will  also  be  studied.  Optimal  trajectories  for  the  case  where  the 
final  coast  to  the  target  is  replaced  by  a  thrusting  stage  will  also  be  determined.  The  problem  of 
optimal  intercept  from  a  space  base  will  be  extended  to  the  case  where  the  carrier  is  an  ascending 
rocket,  namely,  the  (n-l)th  stage,  and  to  the  case  where  the  carrier  is  a  supersonic  aircraft 
maneuvering  at  high  altitude. 

Additional  research  on  methods  of  guidance  and  control  of  the  final  interceptor  stage  will 
consider  the  use  of  more  accurate  target  and  interceptor  prediction  methods,  the  effect  of  guidance 
measurement  errors,  and  the  effect  of  uncertainty  in  target  thrust  termination  time.  The  influence  of 
each  of  these  effects  on  the  performance  attainable  with  the  single-jet  attitude  control  scheme  will 
be  studied. 
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A  great  deal  of  work  on  hypersonic  flow  was  done  in  the  period  1955-1970,  with  some 
attention  paid  to  boundary  layer  blowing.  Although  the  greatest  part  of  this  work  had  to  do  with 
the  flow  over  blunt  bodies,  there  were  some  studies  involving  slender  bodies.  More  importantly, 
multi-component  flows,  dissociation  of  the  air,  gas  reactions  involving  effluents  from  the  surface, 
and  the  strong  effects  of  heat  transfer  received  much  attention  in  those  analyses.  Because 
numerical  techniques  were  in  their  infancy,  analytical  methods  and  the  necessary  attendant 
simplifying  assumptions  were  used.  In  the  meantime,  the  use  of  computational  methods  has  come 
to  the  forefront  and  more  realistic  closure  models  have  allowed  the  numerical  analysis  of  turbulent 
boundary  layers.  Nevertheless,  the  literature  has  much  to  offer  for  this  problem  and  so  the  First 
few  months  of  this  effort  were  used  in  a  literature  search.  Particular  emphasis  was  given  to  papers 
on  gas  injection  into  boundary  layers  in  hypersonic  flows,  or  at  least  compressible  flows,  on 
hypersonic  flow  over  slender  bodies,  and  particularly  over  wedges,  and  on  multi -component  flow 
fields,  and  in  particular,  the  calculation  of  multi-component  transport  properties. 

The  problem  considered  is  that  of  two-dimensional  hypersonic  flow  over  a  slender  wedge, 
from  the  upper  surface  of  which  a  gas,  not  necessarily  air,  is  blown.  Figure  2.1  contains  sketches 
of  the  case  where  the  blowing  is  distributed  over  the  whole  surface  (2.1a)  and  the  case  where  it  is 
confined  to  strips  (2.1b),  and  illustrates  the  notation  employed.  It  should  be  noted  that  s(x)  is 
indicated  to  be  the  distribution  of  effective  body  shape.  Since  it  is  measured  relative  to  the  surface 
of  the  wedge,  it  is  a  displacement  thickness.  For  strong  blowing,  that  is  when  the  viscous 
boundary  layer  is  blown  off  the  wall,  and  moreover,  for  the  case  where  this  viscous  layer  is  thin 
compared  to  the  blown  layer  thickness,  the  displacement  thickness  corresponds  to  the  separating 
streamline;  i.e.,  all  the  blown  gas  flow  lies  within  some  bounding  streamline.  For  other 
conditions,  when  the  viscous  layer  is  not  thin,  a  displacement  thickness  must  be  calculated.  In 
Figure  2.1b,  the  shapes  s i (x>,  S2(x),...  indicate  the  displacement  thickness  distributions  for  the 
flow  resulting  from  all  the  strips  except  the  first  (si(x)),  all  except  the  first  two  (S2(x))  etc.  Just  as 
in  typical  boundary  layer  theory,  one  can  consider  the  thickness  of  a  concentration  layer  consisting 
of  the  region  into  which  the  injected  gas  is  blown  and/or  diffuses.  For  example,  following  usual 
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practice,  the  edge  of  the  layer  might  be  defined  as  the  line  along  which  the  concentration  of  injected 
gas  is  a  very  small  constant;  of  course  this  line  is  not  a  streamline.  The  order  of  magnitude  of  the 
ratio  of  the  thickness  of  this  layer  to  the  displacement  thickness  depends  on  the  case  under 
consideration.  Thus,  for  the  case  where  the  injection  velocity  is  of  the  same  order  as  the  v  velocity 
component  in  the  boundary  layer  and  the  injected  gas  mixes  with  the  boundary  layer  air  (weak 
blowing),  the  order  of  this  ratio  depends  upon  the  Schmidt  number.  Indeed,  *he  effect  upon  this 
ratio  of  the  ratio  of  the  molecular  weight  of  the  blown  gas  to  the  molecular  weight  of  air  may  prove 
to  be  an  important  parameter  in  the  problem  of  obtaining  a  given  pressure  distribution  for  a 
minimum  amount  of  blown  gas.  In  any  event,  the  fundamental  problem  considered  is  that  of 
ascertaining  the  effective  body  shape  s(x)  for  a  given  blowing  distribution  vw(x)  so  that  the 
corresponding  pressure  distribution  on  the  wall  or  body  surface  can  be  calculated;  finally  it  is 
desired  to  solve  the  inverse  problem  so  that  the  blowing  distribution  vw(x)  can  be  found  for  a 
desired  pressure  distribution  pw(x). 

In  hypersonic  flow  without  blowing,  the  shock  wave  on  a  slender  wedge  lies  close  to  the 

body  and  thus  to  the  boundary  layer  over  the  body.  For  some  conditions,  the  inviscid  shock  layer 

downstream  of  the  shock  wave  has  a  thickness  large  compared  to  the  boundary  layer  thickness  and 

for  others  .nay  be  of  the  same  order  of  magnitude,  and  in  fact,  may  merge  with  the  boundary  layer. 

In  addition,  if  the  nose  is  slightly  rounded,  the  shock  wave  is  detached  and  a  thin  high  entropy 

layer  is  found  adjoining  the  wall.  The  parameters  which  control  the  type  of  flow  field  found  are 
the  Mach  number  of  the  free  stream  M^,  the  Reynolds  number  based  on  the  free  stream  properties 

and  the  body  length  Re^,  and  the  wedge  angle  9  and  nose  radius.  Tire  hypersonic  similarity 

parameter,  written  as  in  terms  of  the  parameters  already  defined,  is  taken  to  be  0(1)  or  large 

compared  to  one  for  hypersonic  flow.  The  parameter  which  characterizes  the  relative  order  of  the 
thicknesses  mentioned  above  is  the  viscous  interaction  parameter 


cM3 

ao 


(2.1) 


where  c  is  a  constant.  Thus  for  xi  «  1,  there  is  a  weak  interaction  between  the  inviscid  and 

boundary  layer  flow;  in  other  words,  the  boundary  layer  is  thin  compared  to  the  shock  layer.  As 
XL  increases,  the  boundary  layer  becomes  thicker  relative  to  the  shock  layer  until  a  strong 

interaction  occurs;  then,  they  are  the  same  order  of  magnitude  and  the  solutions  for  each  layer  must 
be  found  simultaneously.  Finally,  for 
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(2.2) 


Xl=0(M 2j 

and  for  xT  »  M  2,  the  layers  are  merged;  they  are  indistinguishable. 

When  gas  injection  is  introduced,  as  an  additional  mechanism  to  change  the  displacement 
thickness,  it  is  still  instructive  to  characterize  the  flow  using  x^-  However,  now  the  ratio  of  the 

order  of  the  displacement  thickness  without  blowing  to  that  with  blowing  S’^/Sb'*'  may  be  used  in 
addition,  as  indicated  in  Figure  2.2,  where  various  physical  problems  and  flow  structures  are 
associated  with  parametric  regimes.  It  should  be  noted  that  there  are  no  definite  limits  on  any  of 
the  regimes.  The  dashed  lines  are  drawn  arbitrarily  simply  to  illustrate  that  there  are  different 
regions  in  the  sketch;  in  actuality,  there  is  a  gradual  change  from  one  flow  regime  to  the  next.  It  is 
instructive  to  picture  the  flow  fields  associated  with  various  parametric  regimes  shown  in  Figure 
2.2.  Several  such  pictures  for  a  flow  over  a  half  wedge  are  shown  in  Figure  2.3,  where  the 
identifying  letters  correspond  to  those  in  the  circles  in  the  sketch  in  Figure  2.2.  For  both  weak  and 
strong  blowing,  as  XL  increases,  merged  layers  (not  shown  in  Figure  2.3)  occur.  As  seen  in 
Figure  2.3,  the  most  striking  feature  of  the  flow  field  with  blowing  is  that  as  the  blowing  velocity 
increases,  the  boundary  layer  can  be  completely  blown  off  the  wall.  It  is  this  feature  which  will  be 
used  to  characterize  the  terms  weak  and  strong  blowing.  Thus,  as  long  as  the  boundary  layer  is 
attached,  the  blowing  is  w-eak,  while  for  strong  blowing  the  boundary  layer  is  blown  off  the  wall. 


It  was  decided  to  attack  the  problem  from  two  different  viewpoints.  First,  one  of  the 
important  problems  in  developing  a  computational  method  of  solution  for  compressible  flow  is  the 
location  and  crispness  of  the  shock  wave.  This  is  a  more  difficult  problem  in  hypersonic  flow  over 
a  slender  body  because  the  shock  wave  is  at  a  relatively  small  angle  relative  to  the  direction  of  the 
undisturbed  flow.  Hence,  work  is  being  done  on  developing  a  code  which  will  handle  flow  over 
slender  bodies  with  arbitrary  shape,  corresponding  to  a  wedge  with  blowing.  Since  the  inviscid 
flow  over  an  effective  body  (physical  body  plus  displacement  thickness)  is  being  considered,  the 
F.uler  equations  are  being  solved.  The  other  point  of  view  under  consideration  has  to  do  with  the 
analyses  of  the  near  wall  layers  when  blowing  takes  place.  It  was  decided  to  consider  the  hard 
blowing  case  first  for  two  reasons.  First,  the  weak  blowing  case  is  essentially  contained  in  the 
hard  blowing  formulation;  i.e.  no  new  mechanisms  need  be  considered.  Second,  one  of  the 
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important  sub-cases  in  the  hard  blowing  formulation  is  that  where  the  viscous  shear  layer,  which 
results  when  the  boundary  layer  is  blown  off  the  wall,  is  thin  compared  to  both  the  blown  layer 
and  the  shock  layer.  Hence  the  problems  is  as  sketched  in  case  (c)  in  Figure  2.3,  and  it  is  seen  that 
for  distributed  blowing  over  the  whole  surface,  the  flow  is  inviscid.  Thus,  solutions  found 
analytically  could  also  be  calculated  numerically,  because  in  this  special  case  the  whole  flow  Field  is 
inviscid;  such  test  problems  can  be  used  to  validate  the  numerical  code  as  well  as  give  useful 
information.  Moreover,  the  extension  to  the  case  where  the  viscous  shear  layer  thickness  is  no 
longer  negligible  is  relatively  simple.  The  work  discussed  in  this  report,  then,  is  concerned  with 
numerical  solutions  to  the  inviscid  flow  equations  (Euler  equations)  and  to  solutions  relating  the 
equivalent  body  shape  to  the  distribution  of  blowing  velocity  at  the  wall  for  the  case  of  strong 
blowing;  in  the  latter  case  problems  with  flow  fields  similar  to  those  shown  in  both  (c)  and  (d)  of 
Figure  2.3  are  considered. 

2. 1  Analytical  Approach 

In  the  analyses  which  follow,  the  lengths  are  made  dimensionless  with  respect  to  a  length  L 
(overbars  denote  dimensional  quantities)  which  for  the  present  is  arbitrary,  but  which  will  be  used 
later  to  denote  the  length  of  the  body.  Moreover,  the  dimension  normal  to  the  flow  is  ordered  by 
5l*.  the  basic  small  parameter  of  the  problem,  and  defined  as  the  value  of  the  nondimensional 
displacement  thickness  at  x  =  1  (x  »  L).  Thus  y  and  the  displacement  thickness  are  denoted  by 


y  =  y/&L  (2.3a) 

yd  =  5l  s(x)  (2.3b) 

where,  as  previously  indicated,  s(x)  gives  the  distribution  of  displacement  thickness.  For 
«  1,  y,i  is  the  dividing  streamline  between  the  blown  flow  and  the  shock  layer  flow;  the 

viscous  shear  layer  of  negligible  thickness  lies  along  yj  also. 

If  P^,  p^,  and  U„  are  the  undisturbed  pressure,  density,  and  velocity  respectively,  then 
asymptotic  expansions  for  the  velocity  components  and  the  pressure  are,  for  the  case  where 
blowing  is  small  enough  that  the  blown  layer  flow  is  incompressible. 


q*=^-(5i)1/2  u(x>  y)  +- 

(2.4a) 

„  _  *  \n  ~ 

qy  =  (\)  v<x’  y)  +- 

(2.4b) 

-  _  _«2*~ 

p-Poa  =  Po.UM  SL  p(x,  y)  +... 

(2.4c) 

Here  qx  and  qy  are  the  velocity  components  in  the  x  and  y  directions  respectively.  If  the  blowing 
velocity  is  large  enough  that  the  blown  layer  flow  is  compressible  then. 


.PE  T  .1/2 

q  =  (— = — ~)  u(x,  y)  +...  (2.5a) 
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q  =  5lv(x,  y)  +...  (2.5b) 

y  m 

W 

p  *  p-  U2  (8.)  p(x,  y)  +...  (2.5c) 

where  jp«,  «  p  and  is  neglected,  R0  is  the  universal  gas  constant  and  Tw  and  mw  are  the  gas 
temperature  at  the  wall  and  the  molecular  weight  of  the  gas  injected  at  the  wall,  respectively. 

At  hypersonic  speeds  the  boundary  layer  along  a  solid  surface  can  have  a  significant 
displacement  effect  on  the  external  inviscid  flow,  an  effect  which  can  be  greatly  increased  by 
surface  blowing.  Blowing  at  any  flight  speed  can  be  strong  enough  that  the  boundary  layer  leaves 
the  surface  and  the  thickness  of  the  "blown  layer"  is  large  in  comparison  with  that  of  the  separated 
free  shear  layer.  This  "blowhard"  problem  has  been  studied  for  several  cases  by  Cole  and  Arocsty 
( 1  ],  using  a  systematic  asymptotic  approach.  Their  work  has  been  chosen  as  the  starting  point  for 
the  analytical  pan  of  the  present  study. 


Cole  and  Arocsty  consider  thin  shapes  and  'how  that  the  pressure  in  the  blown  layer 
depends,  in  a  first  approximation,  only  on  the  streamwise  coordinate  x.  In  most  of  their  examples, 
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the  results  express  the  shape  of  the  dividing  streamline,  which  separates  die  blown  gas  from  the  air 
which  has  passed  through  the  shock  wave  from  the  leading  edge,  in  terms  of  the  blowing  velocity 
and  the  pressure  distribution. 

To  obtain  a  general  idea  of  the  magnitudes  of  the  aerodynamic  forces  available  through 
blowing,  the  first  rather  simple  calculation  concerned  hypersonic  flow  past  a  wedge,  with  a  thin 
ogival  shape  chosen  for  the  dividing  streamline;  i.e.  for  this  calculation,  s(x)  was  chosen,  and  the 
required  vw(x)  and  the  resulting  p(x)  were  calculated..  For  a  blown-layer  thickness  small  in 
comparison  with  the  wedge  thickness,  the  small  pressure  perturbation  is  linear  in  the  slope  of  the 
dividing  streamline,  and  for  constant  surface  temperature  the  density  is  nearly  constant  in  the 
blown  layer.  The  streamwise  velocity  component  in  the  blown  layer  is  then  found  in  terms  of  the 
pressure  from  the  incompressible  form  of  the  Bernoulli  equation.  Integration  of  the  definition  of 
the  stream  function  leads  to  an  expression  for  the  scaled  dividing-streamline  shape  s(x)  in  terms  of 
the  surface  blowing  velocity  vw(x)  and  the  pressure  p(x): 


,  Pw  x1'2 
s(x)  =  ( - ) 

2p 

'  CM 
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{  p(^)  -  p(x)}1/2 


(2.6) 


Inversion  of  the  integral  equation  and  substitution  of  the  linearized  pressure  formula  for  supersonic 
flow  provides  an  explicit  expression  for  vw  in  terms  of  p. 


For  an  ogival  s(x)  it  is  easy  to  evaluate  the  resulting  integral  analytically.  The  results  can 
then  be  concerted  to  a  plot  of  mass  flow  against  force,  as  shown  in  Figure  2.4  for  particular  values 
of  the  parameters.  If  there  were  no  blowing,  the  pressure  force  on  the  wedge  surface  (normalized 
in  the  same  way  as  in  the  figure)  would  be,  in  the  appropriate  limiting  case, 
{(y+-  l)02/2)  x  1(}3  =  23.  The  values  shown  in  the  figure  are  as  large  as  one  half  this  value. 


If  the  blowing  velocity  is  somewhat  larger,  the  blown-layer  thickness  is  no  longer  small  in 
comparison  with  the  wedge  thickness.  The  relative  pressure  change  caused  by  the  blowing  also  is 
no  longer  small,  and  density  changes  in  the  blown  layer  can  not  be  neglected,  so  that  the 
compressible  Bernoulli  equation  is  required  The  scaled  dividing-streamline  shape  is  now 
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(2.7) 


2Y  o  p(x)  p(^) 

where  Hc  =  total  enthalpy.  A  curve  similar  to  that  in  Figure  2.4  is  obtained  by  numerical 
integration,  and  shows  (Figure  2.4)  the  increased  force  available  with  these  larger  mass-flow  rates. 

Additional  curves  for  th.  latter  case  of  a  compressible  blown  layer  are  shown  in  Figure  2.5 
for  a  dividing  streamline  s(x)  =  (const.)  {1  -  (1  -  x)3},  0  <  x  <  1,  and  s’(x)  =  0,  x  >  1.  That  is,  the 
dividing  streamline  becomes  parallel  to  the  wedge  surface  at  x  =  1,  with  continuous  curvature.  It  is 
found  that  the  blowing  velocity  then  approaches  zero  at  x  =  1,  as  seen  in  the  Figure.  If  instead,  as 
in  Figure  2.6,  s(x)  =  (const.)  {1  -  (1  -  x)2},  0  <  x  <  1,  and  s’(x)  =  0,  x  >  1,  the  blowing  velocity 
drops  to  zero  discontinuously.  These  results  are  consistent  with  series  expansions  of  the  integral 
equation  about  x  =  1  which  have  been  carried  out  for  the  two  cases.  These  cases  are  helpful  not 
only  in  illustrating  levels  of  available  forces,  but  also  in  indicating  the  constraints  on  vw  when 
blowing  is  stopped  at  some  point,  as  in  strip  blowing. 

By  contrast,  the  pressure  forces  for  very  weak  blowing  may  be  too  small  to  be  useful.  For 
a  laminar  boundary  layer,  the  displacement  thickness  for  high  Mach  number  is 
5*  »  0  1  Mm2  x/VRex,  where  Rex  is  the  Reynolds  number  based  on  free-stream  quantities  and  x, 

and  a  linear  viscosity- temperature  law  has  been  assumed  for  simplicity.  Without  blowing,  the 
resulting  pressure  perturbation  is  linear  in  the  slope  of  the  equivalent  displacement  surface,  giving  a 
pressure  force  (again  normalized  as  in  Figure  2.4)  of  about  0.3  for  Rex  =  105  and  for  the  wedge 
angle  and  Mach  number  of  Figure  2.4.  With  blowing  this  would  be  multiplied  by  a  numerical 
factor,  say  2  or  3,  and  clearly  is  still  far  smaller  than  typical  values  in  Figure  2.4. 

As  mentioned  above,  some  effort  has  been  spent  on  ascertaining  the  various  sets  of 
conditions  which  can  or  must  exist  as  vw(x)  — >  0  at  a  given  point  on  the  wedge  surface.  Such 
analysis  is  necessary'  in  order  to  handle  properly  the  solutions  for  strip  blowing  and  indeed 
blowing  with  multiple  strips.  It  is  evident  from  physical  considerations  and  illustrated  in  Figures 
2.5  and  2.6  that  over  the  length  of  the  wedge  the  wall  pressure  must  decrease  if  the  gas  in  the 
blown  layer  is  to  go  downstream.  That  is,  the  pressure  gradient  in  the  y  direction  is  negligible  and 
so  the  pressure  disuibution  at  the  wall  is  that  which  holds  throughout  the  blown  layer;  if  the 
pressure  rose,  the  fluid  would  be  accelerated  upstream.  However,  it  may  be  desirable  to  obtain 


pressure  distributions  quite  different  from  those  shown  in  Figures  2.5  and  2.6,  so  various 
distributions  of  vw(x)  must  be  considered.  Strip  blowing  is  a  variation  with  promise  because  it 
appears  to  allow  for  local  increases  in  pressure.  That  is,  from  Figure  2.1b,  with  supersonic  flow 
over  s(x),  it  is  seen  that  at  each  discontinuity  in  s(x)  the  pressure  would  rise  in  a  short  distance, 
then  decrease  again  as  x  increases.  The  result  appears  to  be  several  local  increases  in  p(x)  each 
followed  by  a  decrease;  with  the  average  pressure  decreasing  as  one  goes  from  the  vertex  of  the 
wedge  to  some  downstream  point.  This  will  allow  for  a  quite  different  p  distribution  from  those 
found  from  the  s(x)  and  vw(x)  distributions  given  in  Figures  2.5  and  2.6. 

At  the  end  of  a  strip,  the  manner  in  which  vw(x)  -»  0  affects  the  distribution  of  s(x)  and 
thus  of  p(x).  It  can  be  shown  that  s'(x)  -•  0  there  if  there  are  no  other  mechanisms  to  change  the 
pressure.  From  Equation  2.7,  then,  this  can  be  shown  to  require  that  vw(x)  0  continuously 
rather  than  discontinuously.  The  manner  in  which  vw-»  0  fora  given  flow  field,  i.e.  the  functional 
form  needed  to  assure  that  s’(x)  — »  0  continuously  is  found  using  the  derivative  of  Equation  2.7, 
which  relates  s’(x)  to  vw(x)  and  p(x),  and  another  equation  relating  two  of  the  three  unknowns. 
For  first  calculations,  the  tangent  wedge  approximation  is  being  used.  The  resulting  equation  for 
p(x)  and  the  equation  showing  the  proper  order  of  Qw  are: 

p(x)  =  (•— ^-)  (s'(x)  +  6w  )2  (2.8a) 

ew  =  (8*)9w  (2.8b) 

The  final  relation  with  which  one  works  is  an  integro-differential  equation,  of  compiex  form.  The 
analysis  leading  to  the  functional  form  for  vw(x),  and  thus  for  s(x)  and  p(x),  which  is  physically 
correct  and  gives  non-singular  behavior  at  the  edge  of  the  strip,  as  well  as  the  form  of  the  solution 
at  the  beginning  of  the  next  strip,  is  presently  being  carried  cut. 


The  flows  considered  so  far  are  for  x^  «  1;  i.e.  the  flow  pictures  resemble  that  in  sketch 
(c)  in  Figure  2.3.  As  xL  increases,  the  shear  layer  is  no  longer  of  negligible  thickness,  but  of  ihe 

order  of  the  blown  iayer  and  shock  layer  thicknesses;  i.e.,  all  are  the  same  order  and  a  strong 
interaction  taJces  place.  Then  the  flow  appears  as  shown  in  sketch  (d)  of  Figure  2.3  and  in  more 
detail  at  the  top  of  Figure  2.7.  This  case  .also  occurs  for  the  flow  shown  in  sketch  (c)  of  Figure  2.3 
as  the  vertex  or  k  ading  edge  of  the  wedge  is  approached;  i.e.  the  flow  first  follows  sketch  (d)  for 
x  «  1,  nnd  then  follows  sketch  (c)  for  x  =  0(1).  (n  the  following,  the  case  where  the  shear  layer 
thickness  is  no  longer  of  negligible  thickness  compared  to  the  otb-r  layers  is  considered  for 
hypersonic  flow  over  a  flat  plate,  for  simplicity. 
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For  a  flat  plate,  then,  again  using  a  linear  viscosity-temperature  law,  the  boundary -layer 
thickness  OfM^/VRcx)  is  small  in  comparison  with  the  shock-layer  thickness  0(l/Moo)  when 

x  »  M^6  where  is  a  viscous  length  based  on  free-stream  quantities;  vw  is  the 

kinematic  viscosity  in  the  undisturbed  flow.  In  this  range  the  interaction  of  the  boundary-1  ayer 
with  the  external  flow  is  a  "weak"  interaction,  since  the  boundary- layer  flow  can  be  calculated  first 
and  the  correction  to  the  external  flow  determined  later.  For  x  =  the  boundary 

layer  and  external  flow  must  be  calculated  simultaneously,  because  the  boundary- layer  and  shock- 
layer  thickness  are  of  the  same  order  of  magnitude  in  and  Rex.  in  the  range 

Mj  vco/U00  «  x  «  V00/U00  the  interaction  is  called  a  "strong"  interaction,  and  coupled 

self-similar  solutions  are  available  for  the  boundary  layer  and  shock  layer.  When 
x  =  (XM^  \;00/UIX )  ("merged-layer  regime"),  length  scales  in  both  directions  are  of  the  same 

order,  and  the  different  flow  regions  are  no  longer  distinct. 

As  indicated  above,  for  strong  blowing  it  is  anticipated  that  there  may  be  a  significant  flow 
region  where  the  shear-layer  thickness  can  not  be  neglected.  To  gain  some  understanding  of  the 
accompanying  force  changes,  and  for  later  comparison  with  numerical  calculations,  the  case  of 
strong  interaction  with  strong  blowing  is  being  studied.  The  blown  layer,  shear  layer,  and  shock 
layer  are  distinct,  with  self-similar  solutions  available  in  each  region.  These  can  be  obtained 
separately,  with  a  coupling  arising  primarily  because  the  location  of  one  layer  depends  on  the 
displacement  thickness  of  the  layer(s)  below  it. 


To  illusu  ate  the  nature  of  the  solution,  the  form  of  the  transverse  velocity  component  v  in 
each  of  the  three  regions  is  as  shown  below,  in  terms  <~.f  the  appropriate  similarity  variables.  The 
stream  function  \j/  is  defined  in  the  usual  way  by  dvj7/c)y  =  p  qx;  the  coordinate  x  is  defined  by 
x  =  (0oox/~  J/MJ\  so  that  the  strong- interaction  range  is  Mco*4  «  x  «  1. 


.  U  x  .  1/2  „  _  _  .  .  u  x  .  i 

blov-Ti  layer:  v  =  )  v  { (£),  £  =  M«(”  ;-)/ (—~r  ) 

v  Nf  p  v  M  v  1VT 

OO  M  r  0*1  CM  Od  OO  90 


U  X  .1/2 


.U  X.-t/4=  =  ^ 

shear  layer:  v  --=  ( — -—)  v,(^).  $  =  - ^ 

v  M 


6  V  -  .  .6  ' 


shock  layer:  v  =  ( — ~)  9i(£),  |  =  ( - 4 

v  M6  P  JX  vX 


p  v  M  v  M 

'  90  OO  OO  <U)  i 

U  x 


(2.9a, b) 

(2.9c, d) 
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The  solution  for  the  blown  layer  is  found  from  the  results  of  Cole  and  Aroesty  [lj.  In 
particular,  the  surface  blowing  velocity  vw  =  (const.)  0,*/ Vx  is  found  to  give  surface  pressures 

p  =  (const.)  PooU^/Vx,  consistent  with  the  conventional  strong-interaction  theory  but  with  a 

larger  constant  factor,  which  is  to  be  determined.  The  solution  for  the  shear  layer  is  found  by 
numerical  integration  of  the  compressible  boundary-layer  equations  in  self-similar  form;  velocity 
and  total-enthalpy  profiles  are  shown  in  Figure  2.7.  The  solution  for  the  inviscid  shock  layer 
requires  numerical  integration  of  the  self-similar  hypersonic  small-disturbance  equations,  which  is 
currently  being  carried  out. 

2.2  Numerical  Approach 

The  computation  of  hypersonic  flows,  in  the  Mach-number  range  8-30,  puts  before  the 
computational  fluid  dynamicisf  a  number  of  problems  not  encountered  in  the  lower  Mach-number 
range.  One  problem  is  the  loss  of  accuracy  of  conventional  finite-volume  methods  when 
"capturing"  strongly  oblique  discontinuities;  more  problems  are  encountered  in  marching  to  a 
steady  state,  where  loss  of  positivity  of  certain  state  quantities,  non-uniqueness  of  discrete 
solutions  (recently  found  at  NASA  Langley  Research  Center)  and  general  inefficiency  of  classical 
relaxation  methods  on  vector  computers  all  contribute  tc  slowing  down  or  even  halting  the 
convergence  process.  Many  of  these  problems  have  to  be  addressed  when  computing  the  two- 
dimensional  hypersonic  flow  over  a  wedge  with  surface  blowing,  which  is  the  theme  of  the  present 
work. 


In  the  period  covered  by  this  report,  emphasis  was  put  on  the  question  of  accuracy, 
although  matters  of  computational  economy  are  not  ignored,  as  will  become  evident  below. 


The  representation  of  discontinuities  oblique  to  the  computational  grid  with  high  resolution 
is  a  fundamental  problem  which  at  present  is  being  considered  only  by  a  handful  of  researchers.  It 
would  be  easier  to  ignore  it  and  simply  rely  on  the  capacities  of  today's  supercomputers,  namely, 
by  using  greatly  refined  grids.  Such  a  strategy  has  often  forestalled  advances  in  computational 
algorithms,  especially  in  research  environments  equipped  with  the  latest,  top-of-the-line 
computers,  and  eventually  backfires.  In  the  present  situation  it  is  the  pursuit  of  three-dimensional 
flow  simulations  that  necessitates  the  development  of  high-resolution  algorithms,  since  the  number 
of  nodal  points  per  dimension  drops  significantly. 
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There  are  two  ways  to  improve  the  resolution  of  a  flow  computation  without  unduly  taxing 
the  CPU  budget;  local  grid  refinement,  and  local  reconstruction  of  discontinuities.  These  are 
independent  techniques  that  may  complement  each  other  in  practice,  and  both  are  equally  worthy  of 
research  effort  and  support.  Here,  it  has  been  decided  to  develop  the  latter  technique,  although 
funding  has  been  requested  from  other  sources  for  developing  grid-adaptation  techniques,  in 
parallel  to  this  work. 


A  study  of  the  literature  on  the  subject  of  "jump  recovery"  at  the  start  of  the  period,  and 
discussions  of  the  subject  with  a  few  active  in  the  field  (S.  F.  Davis,  NSWC,  P.  L.  Roe, 
Cranfield,  and  C.  Hirsch,  Brussels)  during  the  summer  have  led  to  the  following  insights: 

(a)  there  are  several  models  of  local  flow  that  can  provide  information  about  strong 
waves  present  in  a  discrete  (finite-volume)  solution; 

(b)  no  one  knows  exactly  how  to  include  this  information  in  a  computational  flow 
algorithm  that  is  stable  and  yields  the  desired  accuracy. 

While  apparently  the  greatest  challenge  is  in  (b),  there  is  still  ample  room  for  ideas  regarding  the 
modelling  of  local  flow  based  on  limited  discrete  data. 

Following  an  approach  previously  indicated  in  Reference  [2],  the  initial  work  was  begun 
with  a  least-squares  analysis  of  the  local  flow  field  based  on  only  two  sets  of  flow  quantities, 
describing  the  averaged  states  in  two  adjacent  cells  of  fluid.  If  these  states  can  be  connected  by  a 
single  oblique  discontinuity,  they  should  satisfy  the  jump  equations 

ViQj  =  [f]  cos  a  +  fg]  sin  a,  (2.10) 

where  Q  is  the  state  vector  (components:  mass,  momentum  and  energy  density),  f  and  g  are  the 
flux  vectors  in  a  cartesian  frame,  a  is  the  direct  angle  of  the  normal  to  the  wave  front,  v  is  the  wave 
velocity,  and  [..)  denotes  a  jump  (see  Figure  2.8).  In  computational  practice,  exact  satisfaction  of 
these  equations  never  happens,  but  a  set  of  values  (a,  V)  may  be  sought  that  minimizes  the  length 
of  the  residual  vector  r. 


r  =  -  V[Q]  +  [fj  cos  a  +  [g]  sin  a 


(3.11) 
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If  the  notation  (...)  is  used  for  an  inner  product,  then 


(r.r)  =  V2([Q],  [Q»  -2V([Q],  [f]  cos  a  +  [g]  sin  a) 

+  ((£]  cos  a  +  [g]  sin  a,  [f]  cos  a  +  [g]  sin  a) 
which,  for  fixed  a,  reaches  a  minimum  for 

V  =  {([Q],  [f])  cos  a  +  ([Q],  [g]>  sin  a}/([Q],[Q]) 
With  this  choice  of  V, 


r  =  {[f]  - 


([Q],  [f] ) 

([Q],  [QJ) 


[Q]}  cos  a+  {[g]  - 


([Q],  [g]) 

([Q].  m 


[Q] }  sin  a 


=  a  cos  a  +  b  sin  a 

Therefore 

min  (r,r)  =  (a,a)  cos2  a  +  2(a,b)  sin  a  cos  a  +  (b,b)  sin2a 
V 

=  y  ((a,a)  -  (b,b))  cos  2a  +  (a,b)  sin  2a  +  j  {(a,a)  + 


Now  \| f  is  defined  such  that 


tan  2y  -  2(a,b)/{(a,a)  -  (b,b) } 


More  specifically, 

sin  2\j/  =  -2(a,b)/h, 

cos  2y  =  -{(a.a)  -  (b,b)}/h, 

h  =  V { (a, a)  -  (b,b)}2  +  4(a,b)2 


Then 


(2.12) 

(2.13) 

(2.14) 

(b,b) }  (2.15) 

(2.16) 

(2.17) 

(2.18) 
(2.19) 
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min  (r,r)  =  -  jcos  2(a  -  \j/)  +  ^  {(a, a)  +  (b,b) }  (2.20) 

V 

which  reaches  its  minimum  value  for 

a  =  \jr  +  K7t  (2.21) 

It  is  convenient  to  adopt 


ct  =  ipr  (2.22) 

as  the  essential  solution;  with  a  =  V  +  k  there  is  only  a  change  of  sign  in  V.  To  find  V,  (2.22)  is 
inserted  into  (2.13). 


There  is  one  degenerate  case,  namely 

[Q]//[fl//[g]  (2.23) 

leading  to  a  =  b  =  0.  This  occurs  for  the  inviscid  flow  equations  when  the  wave  shows  up  only  as 
a  density  (or  entropy)  fluctuation.  Such  a  wave  is  linear,  implying  that  the  wave's  velocity  vector 
is  independent  of  the  fluctuations  it  causes.  In  this  case  the  normal  to  the  wave  front  is  not 
meaningful;  the  flow  angle  remains  the  only  useful  direction.  The  formula  for  a  is  therefore 
modified  such  as  to  yield  the  flow  angle  in  the  degenerate  case  (Eq.  2.23).  Inserting  this  value  into 
Eq.  (2.13)  returns  the  flow  speed. 


Since  entropy  waves  do  not  contribute  to  the  detection  of  a  wave  direction  from  two 
neighboring  states,  one  might  as  well  remove  these  from  the  analysis.  Eliminating  entropy 
variations  not  accompanied  by  pressure  variations  leads  to  a  very  simple  formula  for  a: 


tan  (2a)  =  2[ul  [vj  /( [up  -  [v]2) 


(2.24) 
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Tliis  equation  has  two  solutions: 


tan  oq  =  [v]  /  [u] ,  tan  ct2  =  -[u]  /  [v]  (2.25a, b) 

Here  u  and  v  are  the  x-  and  y-components  of  the  flow  velocity.  The  first  angle  c«i  indicates  the 
direction  normal  to  a  shock  wave,  the  second  angle  a.2  the  normal  to  a  shear  wave.  By  removing 
the  entropy  variations — which  were  irrelevant —  bifurcation  has  been  introduced:  it  is  not  clear 
how  to  switch  smoothly  from  one  case  to  the  other  for  general  data.  In  this  respect  the  analysis  is 
still  incomplete. 


Knowledge  of  ot  enables  a  solution  to  be  found  of  the  problem  of  the  interaction  between 
the  two  states,  called  Ql  and  Qr,  in  the  proper  frame  of  reference,  in  particular,  using  the  proper 
projected  velocities.  The  generic  formula  to  compute  the  flux  of  the  interface  between  two  fluid 
cells  follows  from  the  solution  of  the  Riemann  problem  defined  by  the  different  states.  Solving  the 
Riemann  problem,  however,  requires  an  iterative  procedure,  even  when  the  ideal-gas  law  is 
assumed,  which  explains  the  emergence  of  several  highly  useful  "approximate  Riemann  solvers", 
reviewed  in  [3]  and  [4].  Knowledge  of  V  suggests  a  new  approximate  Riemann  solver,  with  the 
algebra  brought  down  to  the  absolute  minimum: 


f  (Ql.  Qr)  =  \  [  { f(OL)  +  ?(Qr)  }  -  IVI  (Qr  -  Ql)]  (2.26) 

The  sign  a  indicates  that  the  fluxes  are  measured  in  the  direction  of  the  wave  vector.  If  the  states 
Ql  and  Qr  can  be  connected  by  a  single  discontinuity,  this  formula  recovers  the  flux  that  follows 
from  the  exact  Riemann  solution. 


Less  clear  is  how  to  compute  a  flux  in  the  direction  a  +  tt/2,  i.e.  normal  to  the  wave  vector. 
This  flux,  indicated  by  g,  is  needed  for  compounding  the  fluxes  in  an  arbitrary  direction,  in 
particular,  in  the  direction  normal  to  the  interface  separating  cells  L  and  R  (see  Figure  2.9).  The 
simplest  formula  is 


31 


(2.27) 


g(QL,  Or)  =  %  { g(QL)  +  I(Qr)  ) 

but  this  leads  to  central-differencing  of  the  g-fluxes  in  the  final  updating  scheme  and  could  be 
unstable.  Numerical  tests  confirm  this  as  will  be  seen.  Davis  [5]  avoided  this  type  of  instability  by 
averaging  g  over  3  cells  rather  than  2,  as  in  Eq.  (2.27),  according  to  an  algorithm  distinguishing  8 
cases.  Present  experimentation  is  aimed  at  finding  a  simpler,  yet  robust  algorithm  using  a  minimal 
number  of  data. 


The  above  analysis  is  not  restricted  to  updating  schemes  of  a  particular  order  of  accuracy. 
In  fact,  tests  of  the  numerical  resolution  of  discontinuities  arc  best  carried  out  on  the  basis  of  a 
first-order  upwind-differencing  scheme,  since  the  penalty  on  ignoring  the  direction  of  a 
discontinuity  is  greatest.  Results  for  an  oblique  shock  show  the  excessive  smearing  typical  for  a 
standard  upwind  scheme.  In  contrast,  a  shock  aligned  with  the  grid  has  only  1  to  2  cells  across. 


For  the  actual  computations  a  scheme  with  second-order  accuracy  is  needed  (one  of  the  so- 
called  K-schemes  tested  in  [6]);  this  is  presently  being  tested.  The  time-marching  experimentation 
will  be  carried  out  with  new  explicit  methods  developed  in  parallel  by  some  other  doctoral  students 
in  the  department.  These  comoine  a  local  but  matrix-valued  preconditioner  with  a  Runge-Kutia- 
type  updating  scheme  that  hides  a  low-pass  filter.  Such  a  marching  algorithm  is  needed  for  a 
successful  multigrid  strategy,  to  be  added  later.  Both  the  preconditioner  and  the  Runge-Kutta 
scheme  depend  crucially  on  the  availability  of  directional  information  such  as  the  angle  a  derived 
above.  The  choice  of  the  marching  scheme,  viz.  explicit  rather  than  implicit,  is  motivated  by  the 
availability  of  vector  computers,  for  which  such  schemes  are  pre-eminendy  suited. 
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The  numerical  work  was  started  by  developing  a  set  of  grids  on  which  to  carry  out  the 
numerical  computations.  The  logical  first  choice  was  a  cartesian  grid  positioned  along  the  wedge 
face  extending  from  the  nose  back  (Figure  2.10).  However,  once  blowing  is  applied  at  the  nose  it 
is  anticipated  that  the  effective  body  shape  will  become  more  blunt  and  the  shock  may  stand  out  off 
from  the  body  and  thus  fall  off  this  grid.  Moreover,  since  the  nose  area  is  anticipated  to  be  a 
critical  area,  a  local  concentration  of  points  may  prove  worth  while.  With  this  in  mind  a  standard 
polar  grid  and  a  C-tvpe  grid  were  developed  (Figs.  2. 1 1  and  2. 12).  The  C-type  grid  in  this  case  is 
the  result  of  taking  a  regular  C  grid  and  letting  the  radius  of  curvature  go  to  zero  at  the  nose.  This 
creates  a  set  of  triangles  at  the  nose  which  should  not  be  a  problem  with  a  Godonov-type  finite- 
volume  scheme.  The  C-type  grid  (as  is  the  cartesian  grid)  is  created  by  laying  out  points  along  the 
body  and  then  constructing,  row  by  row,  nearly  orthogonal  cells  with  specified  wall  heights  By 
controlling  these  wall  heights,  one  can  expand  or  contract  parts  of  the  grid  at  will  (i.e.  concentrate 
points  in  a  specific  area). 


In  order  to  study  various  numerical  effects,  a  series  of  test  problems  are  derived  from  the 
full  problem.  The  test  problem  considered  first  is  the  resolution  of  a  shock  oue  to  a  wedge  with  a 
10°  half  angle  in  Mach  5  flow.  The  numerical  calculations  are  carried  out  on  a  cartesian  grid 
positioned  along  the  wedge  face.  This  results  in  a  shock  at  an  angle  of  about  9.5°  relative  to  the 
lateral  lines  of  the  grid  (exact  solution).  Initially,  a  first-order  upwind  scheme  was  tested  using  a 
Harten-Lax-Roe  approximate  Riemann  solver  for  each  dimension  separately. 


f(QL.  Or)  =  y  l[f(0L)  +  f(QR)J-ivxi(QR-QL)i  (2.28a) 

g(Qr.  Qb)  =  7  IIg(Qr>  +  g(QBJ]-'Vyi(Qr-QB)}  (2.28b) 

where  Vx  and  Vy  are  based  on  one -dimensional  information: 


Vx  =  <[QJ,  [f])/([Q],  [QD, 


(2.29a) 


Vy  =  ([Q],  [g])/([Qj,  [Q]). 


(2.29b) 


This  application  of  an  essentially  "one  dimensional”  scheme  to  two  dimensions  results  in  excessive 
smearing  of  the  shock,  as  seen  in  Figs  2.13  and  2.14.  Figure  2.13  is  a  contour  plot  of  the  density; 
the  x-axis  lies  along  the  face  of  the  wedge,  the  y-axis  normal  to  it.  The  nose  of  the  wedge  is  at  the 
origin.  Figure  2.14  is  a  3-D  plot  of  the  density  jump.  Again  the  nose  of  the  wedge  is  at  the  origin, 
but  now  the  y-axis  lies  along  the  face  of  the  wedge,  the  x-axis  normal  to  it.  The  x-y  plane  is  the 
computational  grid  and  the  magnimde  of  the  density  is  plotted  normal  to  it  along  the  z  axis. 


Next  the  first-order  code  was  extended  to  calculate  the  cell  interface  fluxes  based  on  fluxes 
measured  normal  to  and  parallel  to  the  wave.  The  flux  normal  to  the  wave  is  described  by  the 
formula 


f(QR,  Ql)  =  2  Uf(QR)  +  f(QrJ]-lVKQR-QL)}  (2.30) 

With  V  obtained  from  Eqs.  (2.13  -  2.22),  the  flux  tangent  to  the  wave  is  found  by  the  simple 
averaging  procedure,  Eq.  (2.27).  This  results  in  an  unstable  scheme.  The  calculation  was  started 
with  an  initial-value  distribution  similar  to  that  of  the  exact  solution.  Figure  2.15  shows  the 
distribution  after  3  time  steps.  Figure  2.16  shows  the  distribution  after  233  time  steps  when  one 
cell  pressure  has  become  negative. 


By  adding  dissipation  along  the  wave,  i.e.  using 

£(Qr.  Ql)  =  2  {lg(QRH  £(Ql)HVI(Qr-Ql)}  (2.31) 

the  scheme  becomes  stable  but  the  lower  quality  of  the  non-rotated  algorithm  is  recovered,  as 
shown  by  the  smeared  shock  in  Fig.  2.17. 


With  this  matter  still  unresolved  a  second-order  upwind  scheme  was  implemented.  This 
produces  sharper  shocks  to  begin  with,  since  it  is  assumed  that  the  distribution  of  state  quantities  in 
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each  cell  is  linear  rather  than  uniform,  as  for  the  first-order  scheme.  The  gradient  of  this 
distribution  in  each  cell  is  obtained  by  central  differencing,  e.g. 


/•3u\  A.,„u  +  4.  ,„u 

(~).  =  - JHH—,  (2.32) 

3x  J  2Ax 

but  this  value  is  "limited"  in  order  to  prevent  numerical  oscillations.  It  turns  out  that  small 
oscillations  nevertheless  remain  present;  these  are  sensitive  to  the  strength  of  the  limiter  as  well  as 
to  the  obliqueness  of  the  shock  with  respect  to  the  grid. 

The  limiter  used  for  the  results  of  Fig.  2.18  is  the  weakest  possible  one,  due  to  Van  Leer 
[7]: 


(i  A .  ...  u+A.1A,u 

—  minrnod  (2A  u,  - _J - ,  2A.  2u) 

Ax 

if  sgn  (A.  I/2  u)  -  sgn  (Aj+1/2  u) 

0  otherwise  (2.33) 

For  Fig.  2  19  the  strongest  possible  limiter  was  used: 


minrnod  (Aj  1/2u,  A.+1/2u) 
if  sgn  (Aj  )/2u)  =  sgn  (A.+1/2u) 


V.  0  otherwise  (2.34) 

still  leaving  some  oscillations.  Adjusting  the  mesh  ratio  Ay/Ax  such  that  the  shock  runs  diagonally 
across  the  grid  (see  Fig.  2.20)  removes  almost  all  oscillations  when  the  limiter  in  Eq.  (2.34)  is 
used.  Further  experimentation  is  needed  to  ensure  that  the  combination  of  flux  formulas  and 
gradient  limiter  gives  a  monotone  shock  profile  regardless  of  the  mesh  ratio. 
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b)  Strip  Blowing.  (Dark  strips  on  wall 

indicate  regions  where  gas  is  injected.) 


FIG  2.1  SKETCHES  OF  DISTRIBUTED  AND  STRIP  BLOWING 
FOR  FLOW  OVER  A  SLENDER  WEDGE. 
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FIG  2.2  CHARCTERiZATICN  OF  FLOWFIELD 
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FIG  2  3  SKETCHES  OF  'T.OWFIELDS 
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COMPARISON  OF  BLOWHARD  SOLUTIONS  FOR 
COMPRESSIBLE  AND  INCOMPRESSIBLE  INNER  LAYERS 
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Fig.  2.4  INJECTED  MASS  FLOW  RATE 
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FIG  2.8  Definition  of  a  and  v 
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oining  cells;  both  t  and  g  are 
npute  the  flux  F  normal  to  the 
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THETR=10  DEG. 
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Fig,  2,10  Cartesian  grid  over  wedge; 
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Fig.  2.12  C  type  grid  over  wedge; 
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Fig.  2.14  Plot  of  the  resulting  density  distribution 
when  the  Har ten-Lax-Roe  approximate  Riemann  solver  is 
applied  to  each  dimension  separately.  The  Y  axis  lies 
along  the  wedge  face,  the  X  axis  lies  along  the  normal 
to  the  wedge  face,  and  the  density  is  plotted  along  the 
2  axis  . 
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Hg.  2.15  Plot  of  the  density  for  the  case  when  numerical 
dissipation  is  added  only  normal  to  the  wave.  The  scheme 
has  been  stopped  after  3  time  steps  to  show  the  onset  of 
instabi 1 i ty . 


Fig,  2.1/  Plot  of  the  density  for  the  case  where  numerical 
dissipation  has  been  added  both  normal  to  and  along  the  wave. 
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F:p.  2.18 

Plot  showing  tne  magnitude  of  the  density  cs  found  by  the  second  order  scheme 
using  the  fksx  limiter  as  in  r.q. 


Fig  2,i9 

Plot  showing  the  magnitude  of  the  density  as  found  by  the  second  order  scheme 
using  the  fiux  limiter  as  in  Eq,  2. 34, 


Fig.  2.20 

Plot  showing  the  magnitude  of  the  density  as  found  by  the  second  order  scheme 

using  the  flux  limiter  as  in  Eq.  2.34,  but  evaluated  on  a  grid  scaled  such  that  the 

shock  passes  through  each  cell  at  a  45°  angle  relative  to  the  cell  walls  M  =59=  10° 

°  CD 
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